function b = assemble_rhs(dofs,V,T,d,f,caseNum,varargin)

nt = size(T,1);
% shapes and shapes_deri
[qw,qp] = quad_rule(13);
% reserve spaces
d_max = max(d);
phi_cell = cell(d_max,1);
phi_u_cell = cell(d_max,1);
phi_v_cell = cell(d_max,1);
% pre-set the value of the cell matrix
for i=1:10
     [phi_cell{i},phi_u_cell{i},phi_v_cell{i}] = eval_bb(i,qp(1,:),qp(2,:));
end
% note that must be evaluate at master element ,so use qp(:,:);
dim = max(max(dofs))-min(min(dofs)) + 1;
b = zeros(dim, 1);
for tri = 1:nt
  degree = d(tri);
  m = (degree+1)*(degree+2)/2 ;
  phi = phi_cell{degree};
  v1=V(T(tri,1),:);v2=V(T(tri,2),:);v3=V(T(tri,3),:);
  x13 = v1(1) - v3(1);  %a
  x23 = v2(1) - v3(1);  %-c
  y23 = v2(2) - v3(2);  %d
  y13 = v1(2) - v3(2);  %-b
  J = x13*y23 - x23*y13;  

  qp_glob = zeros(size(qp));
  qp_glob(1,:) = v3(1) + x13*qp(1,:) + x23*qp(2,:);
  qp_glob(2,:) = v3(2) + y13*qp(1,:) + y23*qp(2,:);  
  % evaluate at integrate points, 
  fvec = feval(f, qp_glob(1,:), qp_glob(2,:),caseNum,varargin{:});

%  then assemble the rhs:
  b_loc = zeros(m,1);    % b_loc = J*(((ones(m,1)*fvec).*shapes(qp))*(qw'));
  for i=1:m
     b_loc(i) = (fvec.*phi(i,:))*(qw')*J;
  end
  loc_dof = get_tri_dof(dofs,tri);
  b(loc_dof) =  b(loc_dof) + b_loc; 
end
